# SI Figure 5
# Range of autocorrelation coefficients
rm(list=ls())
a <- c(0,.2,.4,.6,.7,.8,.85,.9,.95,.96,.97,.98,.99)

# Length of time series
n=200

# Economic assumptions 
gdp = 80000 # initial aggregate consumption (billion US dollars)
g = 0.019 # annual growth rate of undisturbed aggregate consumption
eta = c(0,1.45) # marginal utility of consumption
rho = c(0.04255,0.015) # pure rate of time preference (1.5% in DICE-2016)

# Demographic assumptions
pop = 7.5 # billion people
pop.asympt = 11.5 # Asymptotic population, in billion
param = 0.134/5 # DICE-2016 "Population 2050 parameter"
pop.ts <- rep(NA,n)
pop.ts[1] <- pop
for (i in 2:n) { pop.ts[i] <- pop.ts[i-1]*(pop.asympt/pop.ts[i-1])^param}
# for (i in 2:n) { pop.ts[i] <- pop.ts[i-1]}
pop.ts <- matrix(rep(pop.ts,length(a)),ncol=length(a),nrow=n)

# Undisturbed output per capita
gdp.ts <- gdp*(1+g)^seq(1,n,1)
gdp.ts <- matrix(rep(gdp.ts,length(a)),ncol=length(a),nrow=n)
gdp.ts <- (gdp.ts/pop.ts)/1000 # undisturbed consumption per capita expressed in thousands of dollars.

# Discount factors
df.2 <- 1/((1+rho[1])^seq(0,n-1,1))
df.2 <- matrix(rep(df.2,length(a)),ncol=length(a),nrow=n)

s <- 1

# Calculate expected net present value of consumption for a sample of temperature paths differing only in autocorrelation.
rm("gdp","i","param","pop","pop.asympt")
library(doParallel)
sample = 10000
cl <- makeCluster(2)
registerDoParallel(cl)
figure2a <- foreach(j=1:sample, .combine='rbind') %dopar% {
	y <- matrix(nrow=0,ncol=length(a)+5)
	for (sd in c(.2,.4)) {
		x <- matrix(nrow=0,ncol=length(a)+4)
		for (temp in c(2,4,6,8)) {
			ar <- matrix(NA, ncol=length(a), nrow=n)
			for (k in seq_along(a)) {
				ar[,k] <- as.vector(arima.sim(n=n, list(order=c(1,0,0), ar=a[k]), sd=sd*sqrt(1-a[k]^2), n.start=1)) + temp # Initialise at fixed value
			}
			
			# Output
			x <- rbind(x,c(colSums((((gdp.ts*1000*(1/(1+((ar)/20.46)^2 + ((ar)/6.081)^6))) * df.2 * (pop.ts*(10^9))))),temp,eta[s],rho[s],g))
		}
		y <- rbind(y,cbind(x,sd))
	}
	return(y)
}
stopCluster(cl)
colnames(figure2a) <- c(1:length(a),"Temp","eta","rho","g","sd")
figure2a <- data.frame(figure2a)
figure2a[,1:length(a)] <- sum(gdp.ts[,1]*1000 * df.2[,1] * (pop.ts[,1]*(10^9))) - figure2a[,1:length(a)] # Damages

s <- 2

# # Calculate expected net present value of consumption for a large sample (100,000) of temperature paths for each temperature-variance-autocorrelation-tuple.
# rm("gdp","i","param","pop","pop.asympt")
# library(doParallel)
# sample=1000
# for (l in 1001:1100) {
# 	cl <- makeCluster(2)
# 	registerDoParallel(cl)
# 	figure2 <- foreach(j=1:sample, .combine='+') %dopar% {
# 		y <- matrix(nrow=0,ncol=length(a)+5)
# 		for (sd in c(.2,.4)) {
# 			x <- matrix(nrow=0,ncol=length(a)+4)
# 			for (temp in c(2,4,6,8)) {
# 				ar <- matrix(NA, ncol=length(a), nrow=n)
# 				for (k in seq_along(a)) {
# 					# arn <- rep(NA,n)
# 					# adjusted.wn <- rnorm(n,mean=0,sd=sd*sqrt(1-a[k]^2))
# 					# arn[1] <- temp
# 					# for (i in 2:length(adjusted.wn)) {
# 					# 	arn[i] <- temp*(1-a[k]) + a[k]*arn[i-1] + adjusted.wn[i]
# 					# }
# 					# ar[,k] <- arn
# 					ar[,k] <- as.vector(arima.sim(n=n, list(order=c(1,0,0), ar=a[k]), sd=sd*sqrt(1-a[k]^2), n.start=1)) + temp # Initialise at fixed value
# 					# ar[,k] <- as.vector(arima.sim(n=n, list(order=c(1,0,0), ar=a[k]), sd=sd*sqrt(1-a[k]^2))) + temp # Initialise at random value
# 				}
# 				# Output
# 				x <- rbind(x,c(colSums((((gdp.ts*(1/(1+((ar)/20.46)^2 + ((ar)/6.081)^6)))^(1-eta[s]) - 1)/(1-eta[s])) * df.2 * (pop.ts*(10^9))),temp,eta[s],rho[s],g)) #
# 			}
# 			y <- rbind(y,cbind(x,sd))
# 		}
# 		return(y/sample)
# 	}
# 	stopCluster(cl)
# 	name <- paste("sample", l, sep="")
# 	save("figure2",file=paste("~/Documents/Projects/IAM_Stochastic/SIfigureautocorrelation/", name, ".Rdata", sep=""))
# 	rm("figure2", "cl","name")
# 	gc()
# }
#
files <- list.files(path="~/Documents/Projects/IAM_Stochastic/SIfigureautocorrelation/")
library(doParallel)
cl <- makeCluster(2)
registerDoParallel(cl)
x <- foreach(j=1:length(files), .combine='+') %dopar% {
	load(paste("~/Documents/Projects/IAM_Stochastic/SIfigureautocorrelation/", files[j], sep=""))
	return(figure2/length(files))
}
stopCluster(cl)

figure2 <- x
figure2 <- as.data.frame(figure2, stringsAsFactors=F)
figure2$Damage.function <- c("Weitzman")
figure2[,-ncol(figure2)] <- sapply(figure2[,-ncol(figure2)], as.numeric)
colnames(figure2) <- c(1:length(a),"Temp","eta","rho","g","sd","Damage.function")


## Scaling factors for plotting
UofC <- ((gdp.ts[1,1])^(1-eta[s]) - 1)/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar <- UofC-(((gdp.ts[1,1]-0.001)^(1-eta[s]) - 1)/(1-eta[s])) # Utility value of the marginal dollar
UofC <- sum((((gdp.ts[,1])^(1-eta[s]) - 1)/(1-eta[s]))*(pop.ts[,1]*10^9)*df.2[,1]) # NPV of utility of undisturbed aggregate consumption

# NPV of utility of consumption with damages (UofCwD) with deterministic temperatures
x <- matrix(nrow=0,ncol=5)
for (temp in c(2,4,6,8)) {
	ar <- rep(temp,n)

	# Output
	x <- rbind(x,c(sum((((gdp.ts[,1]*(1/(1+((ar)/20.46)^2 + ((ar)/6.081)^6)))^(1-eta[s]) - 1)/(1-eta[s])) * df.2[,1] * (pop.ts[,1]*(10^9))),temp,eta[s],rho[s],g))
}

UofCwD <- x
UofCwD <- as.data.frame(UofCwD, stringsAsFactors=F)
UofCwD$Damage.function <- c("Weitzman")
UofCwD[,-ncol(UofCwD)] <- sapply(UofCwD[,-ncol(UofCwD)], as.numeric)
colnames(UofCwD) <- c("UofCwD","Temp","eta","rho","g","Damage.function")

# Plot
quartz(width = 7, height = 4, type = "pdf", "SI figure autocorrelation", file = "SI_figure_autocorrelation.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,1,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(3,3), heights=c(3,3))

plot(c(0,1),c(0,2*10^14),xlim=c(0,1),type="n",bty="L",xaxs="i",yaxs="i",ylab="",xlab="",xaxt='n',yaxt='n')
for (i in c(2,4,6,8)) {
	r <- apply(figure2a[figure2a$Temp==i & figure2a$sd==.4,1:length(a)], 2, FUN=sd)
	lines(a,r, lty=1)
}
axis(side=1, at=c(0,.2,.4,.6,.8,1), label=c("0",".2",".4",".6",".8","1.0"), tick = TRUE)
axis(side=2, at=seq(0,2*10^14,.5*10^14), label=seq(0,200,50), tick = TRUE, las=2)
mtext(side=3, "SD of NPV damages", outer=T, line=-2, cex=.9, at=0.05, adj=0)
mtext(side=3, "(Trillion USD)", outer=T, line=-3, cex=.9, at=0.05, adj=0)

mtext(expression(paste("E(", Delta, "T)=2", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = sd(figure2a[figure2a$Temp==2 & figure2a$sd==.4,length(a)]),  adj = NA, padj = NA, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=4", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = sd(figure2a[figure2a$Temp==4 & figure2a$sd==.4,length(a)]),  adj = NA, padj = NA, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=6", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = sd(figure2a[figure2a$Temp==6 & figure2a$sd==.4,length(a)]),  adj = NA, padj = NA, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=8", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = sd(figure2a[figure2a$Temp==8 & figure2a$sd==.4,length(a)]),  adj = NA, padj = NA, cex = .8, col = NA, font = NA)


plot(c(0,1),c(0,2*10^14),xlim=c(0,1),type="n",bty="L",xaxs="i",yaxs="i",ylab="",xlab="",xaxt='n',yaxt='n')
for (temp in seq(2,8,2)) {
	lines(a,((UofC-figure2[figure2$Damage.function=="Weitzman" & round(figure2$Temp,0)==temp & round(figure2$sd,1)==.4, 1:length(a)]) - (UofC-UofCwD[UofCwD$Damage.function=="Weitzman" & UofCwD$Temp==temp,1]))/scalar, lty=1) # Economic loss in dollars
}

axis(side=1, at=c(0,.2,.4,.6,.8,1), label=c("0",".2",".4",".6",".8","1.0"), tick = TRUE)
axis(side=2, at=seq(0,2*10^14,.5*10^14), label=seq(0,200,50), tick = TRUE, las=2)
mtext(side=3, "Risk Premium", outer=T, line=-2, cex=.9, at=0.55,  adj=0)
mtext(side=3, "(Trillion USD)", outer=T, line=-3, cex=.9, at=0.55, adj=0)

mtext(expression(paste("E(", Delta, "T)=2", sep="")), side = 4, las=2, line = -10, outer = FALSE, at = (((UofC-figure2[figure2$Damage.function=="Weitzman" & round(figure2$Temp,0)==2 & round(figure2$sd,1)==.4, 1:length(a)]) - (UofC-UofCwD[UofCwD$Damage.function=="Weitzman" & UofCwD$Temp==2,1]))/scalar)[1],  adj = NA, padj = 0, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=4", sep="")), side = 4, las=2, line = -10, outer = FALSE, at = (((UofC-figure2[figure2$Damage.function=="Weitzman" & round(figure2$Temp,0)==4 & round(figure2$sd,1)==.4, 1:length(a)]) - (UofC-UofCwD[UofCwD$Damage.function=="Weitzman" & UofCwD$Temp==4,1]))/scalar)[1],  adj = NA, padj = 0, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=6", sep="")), side = 4, las=2, line = -10, outer = FALSE, at = (((UofC-figure2[figure2$Damage.function=="Weitzman" & round(figure2$Temp,0)==6 & round(figure2$sd,1)==.4, 1:length(a)]) - (UofC-UofCwD[UofCwD$Damage.function=="Weitzman" & UofCwD$Temp==6,1]))/scalar)[1],  adj = NA, padj = 1, cex = .8, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=8", sep="")), side = 4, las=2, line = -10, outer = FALSE, at = (((UofC-figure2[figure2$Damage.function=="Weitzman" & round(figure2$Temp,0)==8 & round(figure2$sd,1)==.4, 1:length(a)]) - (UofC-UofCwD[UofCwD$Damage.function=="Weitzman" & UofCwD$Temp==8,1]))/scalar)[1],  adj = NA, padj = 0, cex = .8, col = NA, font = NA)

mtext(expression(paste("Autocorrelation of ", Delta, "T", sep="")),side=1,outer=T, line=2.5, at=)
dev.off()
